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We study the equation of state in 2+1 flavor QCD with nonperturbatively improved Wilson quarks 
coupled with the RG-improved Iwasaki glue. We apply the T-integration method to nonperturba- 
tively calculate the equation of state by the fixed-scale approach. With the fixed-scale approach, 
we can purely vary the temperature on a line of constant physics without changing the system size 
and renormalization constants. Unlike the conventional fixed-A''t approach, it is easy to keep scal- 
ing violations small at low temperature in the fixed scale approach. We study 2-1-1 flavor QCD at 
light quark mass corresponding to m,-nlm,p ~ 0.63, while the strange quark mass is chosen around 
the physical point. Although the light quark masses are heavier than the physical values yet, our 
equation of state is roughly consistent with recent results with highly improved staggered quarks at 
large Nt. 

PACS numbers: 12.38.Gc,12.38.Mh 



I. INTRODUCTION 

The QCD equation of state (EOS) at high temperature plays key roles to understand the nature of quark gluon 
plasma (QGP), e.g. as inputs of hydrodynamical description of QGP space-time evolution in heavy-ion collision exper- 
iments [l|. Lattice QCD simulations provides us with the only systematic way to calculate the EOS nonperturbatively 
without resorting to phcnomenological assumptions. 

For a quantitatively reliable evaluation of EOS in QCD, it is indispensable to incorporate dynamical up, down and 
strange quarks. However, dynamical quarks requires much computational cost on the lattice. Most calculations of 
EOS have been made in the fixed-TVt approach, in which the temperature T = {Nta)~^ is varied on a lattice with fixed 
temporal size Nt by varying the lattice spacing a through a variation of coupling parameters on a line of constant 
physics (LCP). Here, we note that a sizable fraction of the total computational cost is required to systematically carry 
out zero-temperature simulations to determine the location of the LCP, to get basic information such as the scale and 
beta functions on the LCP, and to renormalize finite temperature observables such as the EOS at each simulation 
point. In QCD with dynamical quarks, such systematic simulations are quite demanding. 

We adopt the fixed-scale approach, in which we vary T by varying Nt at a fixed a • In this approach, because all 
the simulations are done with the same values of the coupling parameters, they are automatically on the same LCP. 
Furthermore, we need zero-temperature simulation at only one point to renormalize the observables at all T's. Thus, 
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the cost for the zero-temperature simulations can be largely reduced. To take or to confirm the continuum limit, we 
may repeat the calculations at several values of a. As the zero-temperature configurations, we may even borrow high 
statistic configurations on fine lattices, which were generated for spectrum studies at T = and are open to public 
on the International Lattice Data Grid (ILDG) [sj. 

The fixed-scale approach is complemental to the conventional fixed- A^t approach in several respects: In the very 
high T region where T 0{a^^), the fixed-scale approach suffers from lattice artifacts due to the coarseness of the 
lattice in comparison with the typical extent of thermal fluctuations, while in the fixed- A't approach one can 
keep T-^/a = Nt finite even in the high temperature limit. In the fixed-scale approach, the spatial volume of the 
system is kept fixed at all T's, while, in the fixed- iVi approach, the volume becomes quite small at high T's. Large 
spatial volume is important at light quark masses to suppress volume effects in the hadron spectrum and thus in the 
determination of the scale and LCP. At small T's, typically at T < Tpc where Tpc is the pseudo-critical temperature, 
the fixed-scale approach keeps a small a, while the fixed- A^t approach suffers from lattice artifacts due to large a. It 
should be kept in mind here that the fixed-scale approach requires high statistics in the low T region where we have 
a severe cancellation in the observables due to the zero-temperature subtraction procedure at large A^f. Nevertheless, 
we think it worth to take advantage of smaller overall simulation costs with the fixed-scale approach to calculate the 
EOS in 2-1-1 flavor QCD with small discretization errors around Tpc. 

Another point of our study is the choice of the quark action on the lattice. Most lattice studies of hot/dense QCD 
have been done with computationally less expensive staggered-type lattice quarks [3, However, their theoretical 
basis such as locality and universality are not well established. Therefore, to check validity of these results it is 
important to compare the results with those obtained using theoretically sound lattice quarks, such as the Wilson- 
type quarks. See [a-Q for recent studies of QCD thermodynamics with Wilson-type quarks. Systematic study of the 
EOS with Wilson-type quarks has been done so far only in the case of two- flavor QCD [lO, Ell- We extend the study 
to the more realistic case of 2-|-l flavor QCD, using a nonperturbatively improved Wilson quark action coupled to a 
RG-improved Iwasaki gauge action. 

Thanks to the fixed-scale approach, we can take advantage of using the zero-temperature configurations on the 
ILDG. Using the same combination of lattice actions as ours, the CP-PACS-I-JLQCD Collaboration has generated a 
set of zero-temparature configurations in 2-fl flavor QCD and has studied their hadronic spectrum [12, 13). Another 
attractive point of the fixed-scale approach in a study with improved Wilson quarks is that, unlike the case of the 
fixed- A^f approach, we can keep the lattice spacing small at all temperatures and thus can avoid extrapolating the 
non-perturbative clover coefficient csvv to coarse lattices on which the improvement program is not quite justified. 

Choosing a simulation point of the CP-PACS-I-JLQCD Collaboration, we carry out finite temperature simulations 
to perform the first calculation of the EOS in 2-1-1 flavor QCD with improved Wilson quarks. Although the light 
quark masses studied are heavier than their physical values yet, we find that the EOS obtained is roughly consistent 
with recent results using highly improved staggered quarks in the fixed- Aff approach at large values of Nt . 

In the next section, we introduce the T-integration method which enables us to calculate the EOS nonperturbatively 
in the fixed-scale approach. The lattice set-up and the simulation parameters are summarized in Sect. IIIIl Results 
of gauge observables are presented In Sect. IIVI In Sect. |V]the beta-functions are evaluated. Our results on the EOS 
are shown in Sect. IVII and a summary is given in Sect. IVIII Appendix lAl is devoted to a discussion about the choice 
of the interpolation procedure for the T-integration method. Preliminary results of this study have been reported in 

MM- 

II. T-INTEGRATION METHOD 

In conventional studies of EOS in the fixed- At approach, the pressure p is nonperturbatively estimated by the 
"integration method" [l6j : 

y J bo oh y J bo \ ob /sub 

where V is the spatial volume of the system, Z is the partition function, S is the lattice action with the coupling 
parameters b = (/3, Kud, Kg, ■ ■ ■), and (• • • )sub is the thermal average with zero temperature contribution subtracted 
for renormalization. This relation is obtained by differentiating and then integrating the thermodynamic relation 
p — [T /V) In Z in the coupling parameter space of 6. The initial point 6o is chosen in the low temperature phase such 
that p(&o) ~ 0. 

This method is inapplicable in the fixed-scale approach because b is fixed in the simulations. To overcome the 
problem, we have developed the "T-integration method" Q : Using a thermodynamic relation at vanishing chemical 
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FIG. 1: Temperature vs lattice spacing at each Nt- The horizontal dashed line at a ~ 0.07 fm represents the lattice spacing 
in this study. The vertical shaded line represents the approximate location of the pseudo-critical temperature at our quark 
masses. 



potential, 



where e is the energy density, we obtain another non-pcrturbative estimate of the pressure, 

P^TdT'-^, (3) 
Jt„ T5 ' ^ ^ 

with the initial temperature Tq chosen such that p{Tq) w 0. Here, the trace anomaly e — 3p is calculated as 



3p _ _J_ db /dS_\ 



T^V da \db/,. 



(4) 

' sub 



where a{db/da) is a vector of the beta functions which describes the variation of b along the LCP. 

When we vary T along a LCP by varying 6, the integral in ^ is equivalent to that in ([1]) with the integration path 
chosen to be on the same LCP. However, ([3]) allows us to vary T without varying b. In the fixed-scale approach, we 
vary T by varying Nt- Because Nt is discrete, we have to interpolate the data with respect to T to carry out the 
integration of ([3]). The systematic error from the interpolation should be checked. 

In 0, the T- integration method was tested in quenched QCD and was shown that the systematic error from the 
discreteness of T is under control when a is chosen sufficiently small as adopted in spectrum studies. The EOS 
from the fixed-scale approach was shown to be well consistent with that from the fixed- iVt approach with large Nt 
{Nt > 8), except for the high temperature limit where the fixed-scale approach suffers from lattice discretization errors 
as discussed in Sect. HI 
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FIG. 2: Time history of the Polyakov loop measured on finite-temperature lattices. The horizontal axis is the trajectory 
number, which is 0.5 x HMC steps in our simulations. 
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TABLE I: Simulation parameters and gauge observables. The zero-temperature results {Nt = 58) are taken from 13,] by the 
CP-PACS+JLQCD Cohaboration. Temperature T is determined using 1/a = 2.78 GeV (a ~ 0.07fm) [l^. The length of one 
HMC step is 0.5 trajectories for finite-temperature simulations. St is the molecular dynamics time step and "bin-size" is the 
bin size for gauge observables, both in units of trajectories, "traj." is the generated trajectory length {— 2x HMC steps) after 
thermalization of "therm." trajectories, "plaq." and "rect." are plaquette and rectangular loop expectation values. (L) and 
XL are the bare Polyakov loop and its susceptibility, respectively. 



III. LATTICE SETUP 



We adopt a nonperturbatively 0(a)-iniproved Wilson quark action \17] coupled with the RG-improved Iwasaki 
jauge action [3 to simulate 2-1-1 flavor QCD: 



^9 = -/9e|e^o^m^'(^)+e^i^m^'(^)| 



(5) 



-S"? = E ^^xDlyqi, (6) 

f—u,d,s x,y 

= S^,y - Kf E{(1 ~ ^t^^ U^,t.5x+t,,y + (1 + 7m) Ul_f,^^5^^f,^y} - 5^,y CswW) Kf E '^M^^M^ C^) 
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FIG. 3: Polyakov loop expectation value and its susceptibility as functions of T. The left panel shows the bare results, the 
right panel shows renormalized results using the renormalization scheme of [23|. Xl is multiplied by 2 and Xi-ron is multiplied 
by 0.004 to fit into the same scale. Also shown in the right panel is the results of (Lrcn) from the p4 staggered quark action 
obtained at rri^'^ /ttlI'^'"^ = 0.05 in the fixed TVt-approach at TVt = 8 29j, where the horizontal axis is rescaled using ro = 0.5 
fm. 

with K„ = tyd = i^ud- The clover coefficient csw(/3) has been evaluated nonperturbatively by the Schrodinger functional 
method in Hadronic properties have been systematically studied with this action by the CP-PACS, JLQCD and 
PACS-CS Collaborations, down to the physical point (Tl [T9|-[2^ . 

In this study, we use the zero-teroperature configurations by the CP-PACS and JLQCD Collaborations [l^l, which 
are open to public at ILDG/JLDG Q. The CP-PACS-I- JLQCD zero-temperature configurations are available at three 
/?'s, five Kud's and two Ks's, i.e. at totally 30 simulation points. Among them, we choose /? = 2.05, Kud = 0.1356 and 
Kg = 0.1351 which correspond to the smallest lattice spacing and the lightest u and d quark masses {niTr/rnp ~ 0.63) 
with nis near its physical point (m^^^/m^ ~ 0.74). The hadronic radius is estimated to be ro/a = 7.06(3) [23|. Setting 
the lattice scale by ro — 0.5 fm, we estimate the scale as 1/a ~ 2.78 GeV (a ~ 0.07fm). The lattice size is 28^ x 56 
{NgU ~ 2 fm), and the number of thermalized configurations are 650 (6500 trajectories), which are stored every 10 
trajectories. Note that the u and d quark masses are still much larger than their physical values. We are planning to 
extend the study down to the physical point [22,]. 

Adopting the same coupling parameters as the zero-temperature simulation il3], we generate finite-temperature 
configurations on 32'^ x Nt lattices with iVt = 4, 6, • • • , 16. Our generation code is based on the Colombia Physics 
System (CPS) code [13] with the RHMC algorithm for the s quark. We tuned the acceptance rate to be about 80% 
with the HMC step size of 0.5 trajectories. The simulation parameters are summarized in Table ID 

Using the relation between T and Nt, our range of Nt corresponds to the range T = 174-696 MeV at /? = 2.05, as 
shown in Fig. [T] Previous studies of the pseudo-critical temperature Tpc in two-flavor QCD with improved Wilson 
quarks at A^t ~ 6 d, suggest Tpc around 200 MeV for m-^/nip ~ 0.63 in two-flavor QCD. Taking into account the 
effect of the dynamical s-quark and also our larger values of Nt ^ 14 around the pseudo-critical point, we expect a 
smaller value for Tpc. In the succeeding chapters, we show that our data suggests Tpc ^ 190 MeV at our simulation 
point, as shown in Fig. [T]by the vertical shaded line. 

The flxed-scale approach is not applicable at very high temperatures where the lattice spacing a becomes too coarse 
to resolve thermal fluctuations 0. We may estimate a typical length scale of thermal fluctuations by the thermal 
wave length X ^ 1/E where E is an average energy of massless particles at finite T. We then obtain A ~ 1/ (3T) from 
E ~ 3TC(4)/C(3) ~ 2.7r for the Bose-Einstein distribution and E ~ 3rC(4)/C(3) x 7/6 - 3.15r for the Fermi-Dirac 
distribution. Thus, data at T ^ l/{3a) should be taken with care [1^. On the present lattice, the data at T ~ 700 
MeV may suffer from some lattice artifacts. 
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IV. GAUGE OBSERVABLES 

The expectation values of gauge observables are measured every 0.5 trajectories. The results of basic observables 
are summarized in Table HI The time history of the Polyakov loop defined by 

1 1 

^ = yE3Trn^(-.-n.4 (8) 

X T=l 

is shown in Fig. [51 The gauge configurations are stored every 5 trajectories, on which quark observables are measured. 
By examining the bin-size dependence of the errors, we estimate the statistical errors for gauge observables by the 
jackknife method with the bin size listed in TablelH while those for quark observables are estimated with the bin size of 
25 trajectories after thermalization of 1000 trajectories. Static quark potentials measured on the same configurations 
are studied in [2^ [2^. In the followings, we disregard the statistical error in T from that of the lattice scale a, which 
is about 0.5%. Note that, because the scale is common for all T's in the fixed-scale approach, a shift in the scale a 
just causes an overall shift of T. 

The left panel of Fig. [3] shows the results of the Polyakov loop expectation value (L) and its susceptibility xl = 
Nf{{L'^) - (L)2) as functions of T. We find that (L) starts deviating from zero at T - 180-200 MeV, suggesting the 
pseudo-critical point around there. 

For a comparison with the results of previous studies in the fixed- A'^j approach, we have to renormalize (L) . Although 
the additive renormalization constant for free energies are independent of T and thus are common for all T's in the 
fixed-scale approach, the Polyakov loop (L) e~^^^ does receive a T-dependent renormalization. To enable a direct 
comparison with the results of staggered-type quarks, we adopt the renormalization scheme proposed in [27| . i.e. we 
renormalize L such that the singlet free energy from Lren = (-^rcn)^*-^ becomes the Liischer's universal bosonic-string 
potential — 7r/(12r) -I- or at r = 1.5 ro [1^, where a is the string tension at T = 0. Using our potential data at T = 
[23], we obtain Zj-cn = 1.4801(90). Our results for (Lrcn) and corresponding susceptibility Xircn are plotted in the right 
panel of Fig. [31 We note that the dependences on T in these quantities are largely influenced by the renormalization 
factor. In spite of the heavier light quark mass in our study, our results for (iron) agree well with a result from the 
p4 staggered quark action in the fixed- A't approach at TVt = 8 [1^ (see the right panel of Fig. [3]) . Similar agreement 
of (iron) between a smeared Wilson- type quark action and a smeared staggered-type quark action is reported in 0. 

In Fig. [31 we also show the results of Polyakov loop susceptibilities. In the left panel of Fig. [31 we do not see a 
clear peak in xl at 180-200 MeV where Tpc is expected . In the right panel of Fig. [31 a peak of Xircn around these 
temperatures is not excluded, but due to the large errors there. This is in contrast with the case of our previous study 
in quenched QCD adopting the fixed scale approach [s^l, in which we observe a clear peak of xl, and also with the 
cases of full QCD studies adopting the fixed- A'^t approach with staggered-type (see e.g. [3l|) and Wilson- type (gI [25j 
quarks. As a possible cause of the absence of a clear peak in this study, we note that the resolution in T is lower than 
that in our previous quenched study. We may have missed the peak between the simulation points. We also note the 
followings: (i) We probably have a crossover in full QCD around the simulated quark masses instead of the first-order 
deconfining transition in quenched QCD. (ii) Our previous experience with improved Wilson quarks suggests that the 
peak becomes milder with increasing Nt. Our Nt ^ 14 around the crossover point is larger than those adopted in 
previous studies with the fixed- TVt approach, (iii) The aspect ratio Ng/Nt is not large at low temperatures in this 
study. All of these will make the peak milder and thus more difficult to be detected when the present resolution in T 
is not fine enough. 

V. BETA FUNCTIONS 

To evaluate the trace anomaly according to ([4]), we need the beta functions a{d/3/da) and a{dKf/da) (/ — ud 
and s). In this study, we define LCP's by mj^/mp and m^^^/m^ at T = 0. The beta functions are determined 
nonperturbatively through the coupling parameter dependence of zero-temperature observables. We use the data of 
amp, m-^jrap and mr^^^jra^ at 30 simulation points of the CP-PACS-I-JLQCD zero-temperature configurations [l3| 
to extract the beta functions. From a previous experience of two-flavor QCD with improved Wilson quarks in the 
fixed- TV; approach [TT|, we expect that, although a(dK//da)'s are much smaller than a{dj3/da), in the trace anomaly, 
the overall magnitude of the quark contribution proportional to a{dKf / dafs is comparable with that of the gauge 
part proportional to a{dl3/da), but with opposite sign. Therefore, evaluation of the quark contribution is important. 

In our previous attempt [14| , we have tried to evaluate the beta functions by the inverse matrix method, which was 
successful in the case of two-flavor QCD [ll|. In 2-1-1 flavor QCD, we fitted the data of amp, mTr/nip and m^^^/m^ as 
functions of three coupling parameters (/3, K^d, Ks), and inverted the matrix of the slopes of the former in terms of the 
latter to obtain the beta functions. However, it turned out that errors in a{dKf/da) are quite large to calculate the 
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quark part EOS reliably, although the magnitude of the beta functions and the result for the gauge part of the trace 
anomaly are consistent with an expectation from the two-flavor case [l^ . The situation is similar also when we use 
the data of pseudo-scalar decay constant instead of rUp. We find that the large errors in a{dKf /da) are mainly due to 
the matrix inversion procedure, through which all components of the inverse matrix get errors of similar magnitude. 
Because aidnf /da) are much smaller than a{dl3/da), we need more precise values of the slopes to suppress the errors 
in a{dKf/da). In the present case of 2-1-1 flavor QCD, the data points of zero-temperature configurations around the 
simulation point are not dense enough to achieve the required precision of the slopes. 

To avoid the matrix inversion procedure, we now adopt an alternative method, the direct fit method [ll|: We fit 
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the coupling parameters, /3, K„d and as a function of three observables anip^ m,T^jm,p and m^^^/m^. Consulting 
the overall quality of the fits, we adopt the following third order polynomial function of the observables in this study: 
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Note that the fits for three coupling parameters are independent with each other. Figures |4] and [5] show the results 
of the global fit © as functions of 77ipa. The fits with dof — 10 lead to reasonable x^/dof (= 1.63, 1.08, and 1.69 for 
the fit of /?, Kud, and Kg, respectively). 

We define the LCP by fixing m-^/mp and 7(7.^,,/77i0. Then, the beta functions are calculated as ad(3/da = 
[anip) dj3 / d{amp) etc., in terms of the coefficients ci, C2, C5, cs, cio, etc. in The resulting beta functions for 
our LCP {m-t^jmp = 0.6337, rrin^Jm^ — 0.7377) are shown in Fig. [HI as functions of /3. Beta functions at other light 
quark masses are shown in Fig. [T] As the variable to set the scale, we may alternatively adopt ottItt, arriK or amn* 
instead of anip in ([9]). Results of the beta functions, at our simulation point {(3 = 2.05 on our LCP), adopting various 
scale setting variables are listed in Table |IT1 Taking the results from a777p as the central value, we obtain 
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as functions of /3. The simulation point of this study is marked by open square. Horizontal and vertical bars at each data point 
represent statistical errors. 



scale setting 


df3 
da 


xVdof 


dl^ud 

da 


xVdof 


dKs 

da 


xVdof 


amp 


-0.279(24) 


1.6 


0.00123(41) 


1.1 


0.00046(26) 


1.7 


aniTT 


-0.319(21) 


1.2 


0.00179(38) 


0.8 


0.00088(22) 


1.3 


aniK 


-0.252(25) 


1.0 


0.00105(44) 


1.0 


0.00043(32) 


1.3 


aniK* 


-0.215(28) 


1.1 


0.00055(47) 


1.2 


0.00002(36) 


1.8 



TABLE II: Bata functions at our simulation point determined by the global fit ((9| or that with alternative scale setting variables. 
Values of x'^/dof for the fits are also given. 



at our simulation point, where the first brackets are for statistic errors, and the second brackets are for systematic 
errors estimated by the variation of the scale setting. 



VI. EQUATION OF STATE 

With our lattice action ([5]) and ([5]), the trace anomaly (e — 3p)/T^ is given by 

= ) (11) 



T4 7V3 \^ da\dp /^^^ da \ dtiud /,^,b da \ dn, , 
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Nt 


r[MeV] #conf. 




(5'r> 






58 


- 


390 


-4.90487(46) 1.904649(79) 


-4.74878(44) 1.909956(75) 


16 


174 


447 


-0.00380(82) 


-0.00065(13) 


-0.00271(80) 


-0.00052(12) 


14 


199 


447 


-0.0125(10) 


-0.00182(17) 


-0.01007(93) 


-0.00153(16) 


12 


232 


495 


-0.02987(88) 


-0.00443(14) 


-0.02590(86) 


-0.00394(14) 


10 


278 


287 


-0.0448(11) 


-0.00679(18) 


-0.0422(12) 


-0.00646(18) 


8 


348 


319 


-0.0576(11) 


-0.00885(15) 


-0.0592(11) 


-0.00898(15) 


6 


463 


159 


-0.0850(15) 


-0.01394(21) 


-0.0947(15) 


-0.01500(21) 


4 


696 


95 


-0.3216(24) 


-0.04966(38) 


-0.3501(23) 


-0.05266(35) 



TABLE III: Quark contributions to the trace anomaly: = (iVf iVt)"' E,,^ Tr("''''{(l - 7M)(7^,M(i?^)J+/i,:r + (1 + 

'yt^)Ui-i,,i,iD^)~lf,,J and = (A^i^M)"' Ea=.M>. Tr('='''V^,7?.,(L>0^,l- In this table, the zero-temperature results 

{Nt — 58) are raw expectation values, while the finite-temperature results are subtracted by the corresponding zero-temperature 
values. Quark observables are measured every 5 trajectories (10 HMC steps) after thermalization of 1000 trajectories, and their 
errors are estimated adopting the bin-size of 25 trajectories. 

with 




where Nf — 2 ioi f = ud and 1 for f = s. We evaluate the traces in ([12]) and ([T3|) by the random noise method with 
complex U(l) random numbers [l^. The number of noise is 1 for each of the color and spinor indices. Results of the 
quark contributions in (|12p and (jl3p are summarized in Table IIIIl 

In Fig.[8l the results of the trace anomaly (|lip is shown by the solid curve. The curve is drawn by the Akima spline 
interpolation [s^]. The central values are the results using the beta functions with the scale setting variable amp, and 
vertical thin bars represent statistic errors. We repeat the calculation using the values of the beta functions adopting 
alternative scale setting variables to estimate the systematic error due to the beta function. We find that the effect of 
the change of the scale setting variable partially cancels with each other among different beta functions in the trace 
anomaly. Resulting systematic errors are shown by thick vertical bars in Fig. |8l The systematic errors thus estimated 
are smaller than the statistical errors in this study. 

We find that (e — 3p)/T^ is small at T = 174 MeV but shows a peak at T = 199 MeV and decreases towards higher 
T. We note that the peak height of about 7 at T = 199 MeV {Nt = 14) is roughly consistent with recent results of 
highly improved staggered quarks (obtained at Nt = 6-12) in the fixed- -/Vj approach [4, ^5|. The shape of (e — 3p)/T'* 
suggests that Tpc locates between 174 and 199 MeV. The fact that xl shows no peak at T = 199 MeV, as shown in 
Fig. |3l is consistent with the absence of the pseudo-critical point at T = 199 MeV. 

Carrying out the T-integration (jS]) using the Akima spline interpolation for the trace anomaly, we obtain the 
pressure p/T'^ shown in Fig. HI Here, we have chosen the starting point of the integration to be at Nt = 16 where the 
trace anomaly vanishes within the statistical error. The energy density e/T"^ is calculated by p/T'^ and (e - 3p)/T^. 
To our knowledge, this is the first result for EOS in 2-1-1 flavor QCD with dynamical Wilson-type quarks. 

In our previous test in quenched QCD, we confirmed that the choice of the interpolation procedure has only minor 
effects on the EOS [2j] . Because the resolution in T is coarser in the present study, we need to reexamine the influence 
of the interpolation procedures on the final values of the EOS. The results are summarized in Appendix [X] We find 
that the systematic errors due to the choice of the interpolation procedure are small in the EOS in comparison with 
the present statistical errors. 

The overall large errors in p/T'^ and e/T^ are mainly due to the large statistic error in (e — 3p) /T^ atT^ 200 MeV 
— they propagate to higher T's through the numerical integration. The large statistic error in {e — 3p)/T'^ at T ^ 200 
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FIG. 8: Trace anomaly (e — Zp)/T'^, energy density e/T'* and pressure ip/T'^ in 2+1 flavor QCD. The thin and thick vertical 
bars represent statistic and systematic errors, respectively. The curves are drawn by the Akima spline interpolation. 



MeV is caused by the enhancement factor in (fTTj) (note that S is proportional to NtNf). Although the central 
value is largely canceled by the zero-temperature subtraction procedure, the errors are magnified. We find that the 
statistical fluctuation is much larger in the gauge part than in the quark parts. Note that the same difficulty exists 
also in the fixed- iVt approach when we increase Nt towards the continuum limit. In the fixed-scale approach, because 
high statistics is required at very low temperatures only, the overall numerical cost will still be lower than that in the 
fixed- A'^t approach when we try to keep a similar magnitude of discretization errors around the transition temperature. 
In the present test, however, we stop at the current statistics and leave the task for the future investigation at the 
physical point. 

An additional source of errors in Fig. [5] is the spacing of the data points in T: Because our lattice spacing a is 
coarser than that of our previous study in quenched QCD and also because Nt is restricted to be even due to the 
CPS simulation code with the even-odd preconditioning, we cannot have the resolution as achieved in our previous 
study. To improve the resolution in T, we need to develop a simulation code for odd NtS. An alternative way out 
may be to combine results at different lattice spacing a. Note that we can choose small values of a in the fixed-scale 
approach. When a's are well in the scaling region, results for physical observables as functions of T should lie on 
the same curves for these a's, but at different discrete points. After confirming insensibility to a variation of a, we 
may combine the results at different a's to more smoothly interpolate the data in T. We leave application of these 
methods to future studies of EOS at the physical point. 

Besides the large errors, our EOS looks roughly consistent with recent results with highly improved staggered quarks 
near the physical point: The peak of the trace anomaly from the stout quarks locates at T « 190-200 MeV with the 
peak height height of about 4.0 [1]. A preliminary result from the HISQ quarks gives the peak height of about 5.6 at 
T w 200-220 MeV [5]. We recall that our light quark masses are heavier than their physical values. The experience 
with improved staggered quarks suggests that the peak becomes higher as the light quark masses are increased (see, 
e.g., M). 



VII. SUMMARY 

We calculated the EOS in 2-1-1 flavor QCD with improved Wilson quarks adopting the fixed-scale approach 0, with 
which we vary T without varying the system volume on a fine lattice. As the first step towards the EOS with Wilson- 
type quarks in 2-|-l flavor QCD, we made simulations at ra-j^/nip ~ 0.63, taking advantage of the fixed-scale approach 



12 




to make use of high-precision configurations by the CP-PACS+JLQCD Collaboration at T = Ts']. Although the 
light quark masses are heavier than their physical values yet, our EOS looks roughly consistent with recent results 
with highly improved staggered quarks near the physical point 0, Q . 

To extend the study towards the physical point, however, we found a couple of issues to be solved: To obtain 
statistically accurate EOS at low temperatures, we need a large statistics proportional to Nj (a power of Nt is 
reduced due to the average over the lattice sites). This is, however, an unavoidable step to calculate observables 
suppressing discretization errors. Another source of systematic errors in EOS is the limited resolution in T due to 
the discrete variation of Nt in the fixed-scale approach. In the present study, because the lattice spacing a is coarser 
than our previous quenched study, and because Nt is limited to be even due to the simulation program set we have 
adopted, this seems to be non-negligible. To improve the resolution in T, we need simulations at odd values of Nt 
and a finer lattice spacing a. An alternative way will be to combine results at different a's, since we can choose a's 
fine with the fixed-scale approach and thus, after confirming that the discretization effects are sufficiently small in 
the observables under study, we may combine the results at different a's to more smoothly interpolate in T . We leave 
these trials to a forthcoming study with much lighter quarks, adopting the on-the-physical-point configurations by 
the PACS-CS Collaboration [23. 
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Appendix A: Comparison of interpolation procedures 

To carry out the T- integration given by ([3]) , we need to interpolate the data of the trace anomaly at discrete values 
of T corresponding to the discrete values of Nt. In this appendix, we examine the interpolation procedures and their 
influences on the EOS with our data. 

In the left panel of Fig. |9l we apply three different interpolation procedures to our data of the trace anomaly. 
Beta functions with the scale setting variable anip are adopted. Long-dotted line, dotted line and solid line represent 
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the results of straight hne, cubic spline and Akima spline [32*] interpolations, respectively. In our previous study in 
quenched QCD, we have adopted the cubic spline interpolation [2J. With our present data, however, we find the 
oscillatory interpolation curve by the cubic spline interpolation. This is due to the coarseness of the present data 
points — data are available only at even values of Nt. Cubic spline is not stable for data sets with sharp variations. 

In such cases, the Akima spline interpolation [s^ is widely adopted. The Akima spline is a combination of local 
cubic polynomials and is known to suppresses such oscillatory behavior around sharp variations. From Fig. |9l we find 
that the Akima spline leads to a more natural curve smoothly following the data points. Therefore, we adopt the 
Akima spline interpolation in this study. 

To estimate the systematic error due to the choice of the interpolation procedure in the EOS, we perform the 
T-integration with these interpolations. The results for the pressure is shown in the right panel of Fig. [3] The strong 
oscillation of the interpolation curve from the cubic spline is averaged over through the integration, and the results 
of p/T'^ are well consistent for all the three interpolations. We thus conclude that the systematic error in the EOS 
due to the choice of the interpolation procedure is much smaller than the statistical errors. 
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